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An Efficient Algorithm for Spectral Analysis of Heart 
Rate Variability 

RONALD D. BEROER. SOLANOE AKSELROD. 
DAVID GORDON, and RICHARD J- COHEN 

XWM( -W. pr«*« » «<"P* efefcot algorithm f.r the derivation 
of » heart rat, signal from the electrocardiogram. We demonstrate that 
the amplitude spectrum of this heart rate signal mart closely ™.ch« 
hit of the inp^tfe-al «o an Integral pulse frequency modulo »on 
(IPFM model of the heort'l pacemaker than do the spectra of otter 
Serfeed heart ~,e signals. The applicability of this algonthm 
ero»-jpeetral analysis between heart rate and othtr phys.olog.e sign.U 
is ||M discussed. 

• ■ • 1. Introduction 
In a recent paper [I). DeBoer e/ of. compared rwo methods that 
Jp,o\ ^TOy* for the study of heart s vanab. h* In 
a skond^per [2]. the same authors presented an evaluat.no of 
Uwse twomemods and of a third by testing each on a sequence of 
STed RR intervals generated by an integral pulse frequency 
modulation (IPFM) model. An IPFM model .s a devco that in e- 
gr*« its input signal until the result of this integrate reaches a 
l^x threshold, at which point the device sends out a pulse, .resets 
the integtatorto aero, and begins Urn .ntegration anew. Hyndman 
«d M<*n [3J fir^t sugsested the IPFM model as a foocuonal de- 
Srtptfon of tile slno-atrial node, and it remms . useful mode I fee 
* e mechanism by which the autonomic nervous sy stem moduU es 
heart rate- We can represent the operation of an IPf M model math- 
ematically as 

S/» -i 
(1 + m<j»<fr. U> 
n 

Tis the integrator's threshold value, which equals J* * 
each RR interval were there no autonom* modi: a uon of die SA 
node's intrinsic firing rate. The input " L?*& 

where all autonomic influences are lumped together .n *« model 
£. am represented by «M. Obviously. ^^g^ft 
RR interval shonens so that the instantaneous hear, rate vanes tn 
nrooortiontoJtO-r.istheiimeof the ith /t wave. 
T£e?r evaluation of the performance of the various spectral 
Jfi for analysis of he. "«« variability. DeBoer « */. eom- 
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pared the results of each method 10 the spectrum of the i«P««S"»> 
applied to the IPPM model that generated the simulated RR inter- 
vals. They considered the latter as the -true" heart rate spef™™- 
These investigator, demonstrated that all three of the merhods they 
considered introduce significant artifacts that conupt the spectra, 
compared to the ••true" spectrum. ^.siiitv 
For several veats we have been studying heart rare vanab ility 
[4J-(6]. DeBoe'r « al (1) claimed we use one of the 
analys * ihat they presented. We in fact employ an a gomhot that 
E from a1uh.ee methods «ha, DcBocr rial. 12] 
and their caraful analysis of these three methods h«*p^ our 
interest in similarly snalyxing the performance of the bean ra te 
spectral estimation technique we use. In this paper, we present <mr 
computationally efikient algorithm and demonstrate that the spec- 
tral estimate it yields almost exactly matches the "iwe spectrum 
of the input signal to the 1PFM model. 

Kr <t*L labeled the three type* or power spectra! estimates 
that they discussed as 1 ) the spectrum of intervals. 2) the spectrum 
of inverse intervals, and 3) the spectrum of counts. (For a detailed 
description of the methods involved in the computation of these 
spectra, see [M p].) The spectrum of intervals and the spectrum 
of inverse intervals are the discrete Fourier transforms (DFT) 
squared of sequences of numbers corresponding to the RR interval 
durations and their reciprocals, respectively. Each number in these 
sequences corresponds to a single beat; this obviously results tn 
uneven sampling of the process in time. Since these numbers are 
evenly spaced only wnen plotted against beat number, the units of 
the frequency axis of ihese spectra are "cycles per beat" instead 
of "cycles per second." While the number of cycles per beat can 
be converted to an average number of cycles per second by mult.- 
p!yin« the former by to average heart rate, it is not surprising that 
these spectra would differ in appearance from that of the input sig- 
nal to an IPFM model that generated the RR intervals. F*T exam- 
ple, if the IPFM model's input were a sine wave, then there would 
be relatively fewer beats (and thus RR intervals) around the sine 
wave's minima, and relatively more when the sine wave is near its 
rruVrma. Thus, the spectrum of intervals and the spectrum of 
verse intervals are the spectra not of a sampled sine wave, but rather 
of a distorted sinusoid-like signal that appears alternately stretched 
out and compressed. Clearly, such spectra will contain harmonics 
of the fundamental sine wave- 

The speetnjm of counts the Fourier transform squared or a set 
of delta functions on a true lime axis spaced according to the sc- 
ience of RR intervals. This power spectral estimate, denoted 
-*c (/). can be computed analytically as 

Ftif) AT' 1 1 1 

A |*<co*<2rAv>-l) + i sin t2*ft t ) f I G> 

where iV is the number of delta functions in the record and r, de- 
nates the location in time of the *th impulse. (Note that (2) con- 
n'ns terms that compensate for the truncation effects that result 
from computing the Fourier transform of a finite set of delta func- 
tions.) Rompdman rr uf. 17 J have presented a modification ot this 
technique for efficient implementation on a personal computer. 
Since .his spectrum ii that of a true time signal, it is free of har- 
n.ortie artifact* liketho** seen in the spectrum of ,n*rv a l* *nd the 
spectrum of inverse interns. On the other hand, when the inter- 
ns arc generated by tr. IPFM model, the spacing of the delta 
functions (not their amplitude) is modulated by the input signal 
applied to the model. Thus, artifacts will appear in the spectrum 
of counts bi sidebands of the mean repetition rate, as in any fre- 
quency modulated process. 

'lit thi* pap*., all power spectral estimates have Seen normalise by 
lividine by the >qmr* v( ihc J"??" "f !>>c input rignal. 
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Pi. 1 (.1 A seame* of an 6CG signal. tW The hean rare samples^ 
L™lc »^TfiCC signal in <a>. *.«.«.!»«• m f^^, 
^wri&lt RR totatv.ll within the local window centered * r, * «/ 
I andT^ is hlU + e//.. The value of the heart rale at «th ample 
Irfnfu »k« to be ,he number of interval* that fell wiih.n dw local 
Endow centred "thai point divided by <h. width of the window, a. 
d<£rfbJ to *Tte«t- <c> The conesponding instantaneous hear, rate s. B - 
^ThTvU* held during each interval is the rec.proe.1 of the duraUon 
,? ^WvnI.Trhe sample values in (b) are equivalent 
?i.nal tt-V would n*ult from convolution of the signal .n (c) wuh a reef 
angular window thai i* two sample intervals wide. 

11. Description of Algorithm 
In Fia 1 we present a schematic description of the algorithm we 
use to derive a heart rate signal from the ECO signal. The steps 
involved are as follows. First, the ECO is sampled at a sufficiently 
high rate to determine the time locations of the R waves ito what- 
ever accuracy is desired. Ne*t. a sampling rate for the heart rate 
signal is chosen. This is a true frequency (i.e., the heart rate sam- 
ples will be evenly spaced In time ,« this fluency) and may be 
chosen arbitrarily, without regard to the mean heart rate -or the fre- 
quency at which the ECG is sampled. V . . 
defined at each heat rate sample point as the time mterval extendmg 
from the previous sample lo the neat. We then count 
or *R intervals (including fractions thereof) that occur within tha 
local window. Examples of how we compute ™. 
vals are shown in Fig. 1. The value r, of the heart rate a, each, 
sample point is taken to be 

where /, is the sampling frequency of the resulting heart rate signal 
and n,h the number of RR intervals that fell lit the Iocs I window 
centered at the fth sample point. Finally, we estimate the heart rate 
now spectron. from the sequence of heart rate samples. 

Several otaeavatkms may be made regarding thii technique. Ob- 
viously, since the heart rate at e*ch point depends on events M > dte 
ECG both in the recent past and the near Arture, thus method 1 cannot 
be employed in real-time analysis without incurring a -J™ 
heart rale signal produced by nur algorithm may equivalently be 
Wewed « samples or a stepwise continuous instantaneous hart 
S with .rectangular r^f')M^ 

stepwise continuous instantaneous heart rate signal maintains an 
amplitude equal to the reciprocal of the current RR interval for ^e 
duration oMhat** interval [see Fig. 1(c)). Tte»i«rd^nwn 
fradi tanal uchometcr signals (see for example Pig. If>) - UD •<» 
^ he «lue held during the *«h interval. (,.. . ,). « Mt«. . - 
i mi 1/fr. - r. .) As DeBocrcf al. have recently noted |9J. the 
££ J Uomir sign., is flawed on ,wo couni, FM. Ik. 
17gn*na S » the ECG by an entire beat, which may be Inconsequen- 

»Tbe spectrum ean be computed from the sequence of heap, rate samples 
usinTtnv^rs.^n) "imntinn meshnd, |S|. *• J-raJj 

™ntov an FFT-bawd windowed periodngrero melnod: in ihis paper for 
SSI"* resufcofDeBoer^./. 121. we have unload 
the Bartlcn w.ndow. 
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R(f , AmpiiRMje spectra from four different methods used to imlyzc the 
l^t Hgpil to x(r) - 1 + 0.3 cos C2x/Lf>. /„ was 0.15 Ha and the , 
model's threshold ?^ 1.05 s. (a) Spccmmi of ^t^. <>> Spectrum 
to(Mvais. (c) Spectrum of counts, (d) Spectrum of bait rale 

S Intends. The frequency I? C») 0>> w ~ J^*™ f 
units of cycles per sec™! by multiplicadori by the mom ^^^1 
m Scribed in the test- Tike spectium of counts (c) was ««V^«*I 
0?taf» (b), and (d) the dc component was removed, and the amoU- 
mde «Sa were then found by taking the square tool of Aeir cwre- 
^Swcrspeon. which were computed by uSdag ^fastFouner 
Jn^rm ffFTl squared of the time domain sequences. For the heart 
^r^ 

o 0 to and befbK taking the square root Of the <m^P©*tag power 
s£m^!it was first divided by the filter function *T/> 1» (4). 



rial in autospccttzl analysis, but can introduce smfactual phase 
shifts in cross spectra, for example, between heart rate and Wood 
pressure or respiration. Second, the traditional tachometer signal 
provide* a biased estimate of the heart rate since the lowest values 
arc held for too short an amount of rime, and the highest values are 
held for inappropriately tang intervals- The mean heart rate thus 
ar«>ears higher than it should. The nOddeUycd instantaneous heart 
rate signal in Fig, 1(c) avoids both or these complications. 

Convolution of the heart rate signal with the rectangular window 
has the effect on the power spectrum of multiplication by a low- 
pass filter. The shape of the filter HI/) is 



[sin(2T/5ftT a 
L 2*Wr 1 



(4) 



where, again,/, is the Sampling frequency of the 
« thai VU is *e width of the rectangular window in the umo do- 
T& filter passes very little power beyond the ^ 
ft e. f&). and ui effects can be compensated for in the band 0 <: 



f < fJ2 bv mulriplyiQg the power spectrum by \mfh h P»«ice, 
it i£i\ v*%J) inecSTbut consider the speotra) 
£2 only toft < / < /A ^ the - Jf — * W f 
significantly amplifies any aliased power in rhe band//4 < / < // 

Z ' While the choice for/, is arbitrary, we typically chooseA -^ 
to and similarly sample Mood *^£™^ d ^£ 
4 Hz This U an appropriate sampbog rate tor tne stuoy * 

™« Furthermore. *b*a the heart. ra« «njj« «e spaced 
Sto time and are synchronized with thc^plcsof^^ 
^yilogfc signal*. doss^pectnJ e^testawecn «b«e v»n«« 
signals are just as easy to compote f^*"?*^ /» / n for 

It should be noted that a power spectral 
«epw« continuous signal shown in Fig. 1(c) «oJ 
analytically without the oeed to generate samples of the signal usiug 
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Fig. i. Same u Fig. 2. but here the input Stgaal applied to the IPFM model 
wis 1 + 0.3 CO* ilrAO + 0-3 cos (2%fy). where/, = 0.12 Hz and/, 
=t 0.)6 H2. The IPFM model threshold was again 1.05 j. The difference 
in height* of the two peaks at 0.12 ind 0. 16 Hz reflects different degrees 
of anifaerua! power leakage into side lobes of the respective peaks, de- 
spite the presence of equal amounts of power at these two frequencies in 
the input signal. 



the relation 



I 



♦ (cos + 0 - cos (2*A» 
- (sin (2*A-i> " sin j 



(5) 



whsre A 1 is the number of steps in the period of observation and l k 
is the time at the beginning of the Jfcih step. This technique gencr- 
a:es spectral estimates whose only artifacts arc the result of the 
smzll amount or generally high frequency power inevitably present 
tn the discontinuities between adjacent steps. However, h is not a 
practical method because the evaluation of (5) for many differem 
frequencies is computationally very burdensome and cannot be 
made more efficient through the use of FFT-like algorithms. 

Furthermore, efforts to synthesize a discrete heart rate signal 
suitable for analysis with ar. FFT algorithm, by mere sampling 
without low-pass filtering of the instantaneous heart rate signal 
shown in Fig. 1(c), would cause the above-mentioned high fre- 
Cheney artifacts to become aliased into the physiologically impor- 
tam low frequency band of the power spectrum. It is important to 



note thai while our method in effect generates heart rate sample 
values of a piece wise continuous signal thai has undergone the nec- 
essary antialiasing filtering, the particular means by which we 
achieve the filtering operation is very efficient. Our method avoid; 
the computational burden of actual digital convolution. 

IH. Comparison Bktwkkn Qua Methoo anp Other Heart 
Ratk Sphctra 

In order to demonstrate that the spectrum of the heart rate signal 
constructed using our algorithm is relatively free of artifacts, we 
performed the same simulations as DcBocr *r ai. did [X\ by imple- 
menting an IPFM nude! on a digital computer. We then computed 
heart rate spectra first using the three methods they presented and 
then u*ing our*. The spectra shown here are amplitude spectra (i.e., 
the square hkji of the power spectrum) as these accentuate the pres- 
ence of harmonics and other artifacts. 

Pig. 2 shows the results of the simulation where the IPFM model 
input signal was 

rt*l) « 1 + 0.3 cos (lr/ m 0 (6) 
and the modulation frequency f m was 0. 16 Hz and the IPFM thresh- 
old was 1.05 s. 

All four spectra in Fig. 2 show a large peal: at the modulation 
frequency (0. 16 Hz). However, the spectrum of intervals [Fig. 2(a)J 
and the spectrum of inverse intervals [Fig. 2(b)) also contain a sig- 
nificant peak at the first harmonic (0.32 Hi) and a smaller one at 
the second harmonic (0.48 Ha) of the modulation frequency, that 
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arc virtually absent in the heart rale spectrum computed with our 
algorithm [Fig. 2(d)]. Similarly, there is a sideband artifact at 0.472 
Hz (0.952-0.48 Hz) in the spectrum of counts [Fig- 2(c)] that is 
totally absent in the spectrum computed with our algorithm. 

Fig. 3 Shows the results of a second simulation in which the 
input signal applied to the IPFM model was 

j(r) « 1 + 0.3 cos (2*/,/) + 0.3 cos (ixf 2 t). O) 

Here, again, the model's threshold was 1 -05 s. The two modulation 
frequencies /, and/j were 0.12 and 0.16 Hz, respectively. These 
are the same parameters as those used by DeBocr ft af. in their 
second simulation (2]. All of the spectra for this case show large 
peaks at the rwo modulation frequencies. In addition, the spectrum 
of intervals [Fig- 3(a)) and the spectrum of inverse intervals (Fig- 
3(b)] possess artifacts at harmonics of both modulation frequen- 
cies, and the spectrum of counts [Fig. 3(c)) contains sideband ar- 
tifacts at integer multiples Of/i artcl/j away fiu»n the mean rcp»xi- 
tioo rate of 0.952 Hz. Furthermore, the spectrum of inverse 
intervals [Fig. 3(b)] contains a component at 0.04 Hz, the differ- 
ence between me two modulation frequencies. All of these artifacts 
are almost completely absent in the heart rate specirum computed 
using our algorithm [Fig. 3(d)). 

IV. Conclusion 
The preponderance of influences that impinge on heart rate orig- 
inate outside the heart, vary slowly compared to the heart rate, and 
are relatively insensitive to the actual liming of ventricular acti- 
vations. For this reason, we feel that it seems more natural to char- 
acterize heart rate on a real-time axis, rather than against "beat 
number.* ' The IPFM model is consistent with this description since 
it lumps autonomic control and all other factors that affect heart 
rate into a single time-varying signal. Our algorithm provides a 
computationally simple definition of a heart rate signal derived from 
the ECO, and as Figs. 2 and 3 demonstrate, the spectrum of this 
signal very closely matches that of the IPFM mode! input signal. 



In fact, because of the very way we define the heart rate signal, 
were this signal applied as the input to an IPFM device, ihe re- 
sulting sequence of RR intervals would be identical to the sequence 
of RR intervals from which the heart rate signal was derived. 
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